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q Abstract 

We use LVDSMC simulations to calculate the second-order temperature jump co- 
efficient for a dilute gas whose temperature is governed by the Poisson equation with a 
i Qh constant forcing term. Both the hard sphere gas and the BGK model of the Boltzmann 

equation are considered. Our results show that the temperature jump coefficient is dif- 
ferent from the well known linear and steady case where the temperature is governed 
by the homogeneous heat conduction (Laplace) equation. 

00 

1 Introduction 

^ Slip-flow theory is a powerful tool that enables the continued use of the Navier-Stokes de- 

scription as the characteristic flow lengthscale (L) approaches the molecular mean free path 
(A) [14J. It can be rigorously derived from asymptotic solution of the Boltzmann equation 
in the limit Kn = X/L <C 1; such an analysis shows that, in this limit, the Navier-Stokes 
description remains valid in the bulk, but fails near the boundaries [161 E]- Fortunately, 
the kinetic effects associated with the inhomogeneity introduced by the walls are only im- 
portant within a layer of thickness O(A) near the boundaries (known as the Knudsen layer) 
and can be accounted for by a boundary-layer type of analysis where an inner kinetic solu- 
tion is matched to the outer Navier-Stokes solution [TBI [T7] . Slip/jump boundary conditions 
and the associated non- adjustable slip coefficients emerge from this analysis as the matching 
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condition between the inner and outer solution [TBI E7J- Carrying out such an analysis to 
second order in Kn yields second-order slip/jump models [El [FT], which can be very useful 
in a variety of cases |14j . 

Accurate determination of slip coefficients using this rigorous procedure is quite chal- 
lenging in general and becomes increasingly more challenging as the order of the expansion 
increases. Original studies focused on the BGK model of the Boltzmann equation [3J, E], 
for which all first-order and second-order coefficients are known [El H7j- The first-order 
coefficients for the hard-sphere gas have also since been calculated [5]. However, although 
the form of the slip expression is known to second order in Kn, second-order slip coefficients 
for the hard-sphere gas are mostly unknown. 

As the companion paper shows [7J , the recently developed reciprocity relations by Takata 
[U [2] can be used to calculate these coefficients. An alternative approach amounts to ex- 
tracting slip coefficients from hydrodynamic fields by comparing solutions of the Boltzmann 
equation with Navier-Stokes solutions [U [TJ]. In these approaches, in addition to high ac- 
curacy (including low statistical uncertainty if a stochastic method is used for solving the 
Boltzmann equation), care needs to be exercised to avoid comparison of the two solutions 
in the Knudsen layer, where the Navier-Stokes solution is not equivalent to the Boltzmann 
solution [H]. This has led to a number of erroneous results in the past. 

In this paper we use this process to calculate the second-order temperature jump coeffi- 
cient for a dilute gas when the temperature field is governed by the Poisson equation with 
constant forcing term. We calculate this coefficient using the recently developed low- variance 
deviational Monte Carlo simulation method [TJl [El [El HH1 E] , which is naturally suited to 
low-signal problems and thus allows calculations at infinitessimal temperature differences. 
The latter are necessary because finite temperature differences introduce density gradients 
and temperature-dependent transport coefficients which may alter the result. 

Our result is verified and put on a more firm theoretical footing by the companion paper 
[7J which considers a mathematically equivalent time-dependent problem, thus clarifying why 
the temperature jump law and coefficient reported here are in general different from the one 
obtained by linear steady-state analysis [IB] . 

2 Background 

We consider a dilute hard-sphere gas of molecular mass m and molecular diameter a, in 
contact with a planar diffusely reflecting boundary at temperature Tb. We also consider the 
BGK model of such a gas, with collision frequency r _1 . In the case of the hard-sphere gas, 
A = (v^Tmo " 2 ) -1 ; while for the BGK gas A = 2c t/ v /7t, where c = \J2RT is the most 
probable speed based on the reference temperature T , n is a reference number density, 
R = ks/m is the gas constant and ks is Boltzmann's constant. 

The first-order temperature jump condition at the gas-wall interface is given by [El [FT] 
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where T = T/Tq, k = ^-Kn = ^(A/L), |b denotes the boundary location, n is the unit 
(inward) normal direction and L is the characteristic system length scale; the numerical 
constant d\ obtains the non- adjustable values of 2.4001 for a hard sphere gas and 1.30272 for 
a BGK gas [TB]; we emphasize that these values correspond to diffusely- reflecting boundaries. 

The utility of first-order slip/jump models primarily depends on the amount of error 
that can be tolerated. Temperature jump coefficients (both first-order and the second-order 
measured here) turn out to be larger than their velocity slip counterparts. As a result, a 
second-order temperature jump correction becomes important at smaller Knudsen numbers. 
In other words, the first-order result Q is typically adequate for Kn < 0.1. 

Asymptotic expansion to second order in k [TBI H7J for linear and steady problems extends 
([I]) to the following jump condition 



a, dT 
T\ B -T* = d lk - 



+ d 3 k 2 



d 2 T 



dh 2 



(2) 



This condition is valid for a quiescent gas — more precisely, a gas that is quiescent under 
no-slip boundary conditions; in the presence of gas flow, additional terms related to the flow 
stress need to be included [IB] • We also emphasize that according to the analysis that yields 
this condition, for linear and steady problems, energy conservation reduces to 

V 2 f - — — — - (3) 

dx 2 dy 2 dz 2 

where (x, y, z) = (x/L,y/L,z/L), and L is the characteristic problem lengthscale. In the 
special case of one-dimensional problems, equation ^ further reduces to 

d 2 f 

v2T =d^ = °- (4 > 

which makes the value of d% irrelevant. This is actually utilized below to calculate the slip 
coefficient due to a forcing term in the temperature equation. 

In summary, jump condition ^ is to be used when the governing equation is ([3]). Within 
this approximation, d 3 is only known (d 3 = 0) for the special case of the BGK model [I~6l[TT] . 
We also note that Deissler's result [I] for second-order velocity slip and temperature jump is 
based on approximate mean-free-path arguments and does not correspond to a self-consistent 
solution of the Boltzmann equation; as a result, it captures neither the correct form of the 
slip/jump relation nor the correct values of the slip coefficients (e.g. compare equations 
(3.40)-(3.42) in [17] to equations (24a) and (51) in 0). 



3 Calculation of the temperature jump coefficient 

To extract the slip coefficient in a dilute gas governed by the Poisson equation with constant 
forcing term, we simulate the steady state of a one-dimensional gas layer bounded by two 
isothermal, diffuse walls at x = ±L/2 and at temperature To, subject to volumetric heating 
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at a constant rate Q. In dimensionless form, the one-dimensional heat equation with constant 
volumetric heating can be written as 



V 2 T 



d 2 T 
dx 2 



5e 



where x = x/L, 72 is a dimensionless form of the thermal conductivity- 
hard spheres and unity for BGK [T7j — and 



(5) 

equal to 1.9228 for 



LQ 
c Pq 



< 1 



P 



(6) 

noksTo is a 



is the dimensionless form of the volumetric heat addition rate. Here 
reference pressure. 

The asymptotic analysis yielding ^ does not apply to the non-homogeneous equation 
([5]). A rigorous derivation which takes the inhomogeneous term into account by considering 
an equivalent unsteady problem can be found in the companion paper [7j, which shows that 
in a quiescent gas, in one spatial dimension, the resulting second-order slip relation is given 
by 



T 



dT 
dn 



d 2 T 

b on 2 



(7) 



We emphasize that, although the structure of the slip relation is the same as in equation rt2J), 
the second-order coefficient is different. It is also convenient that ^ does not contain (I3; this 
allows calculation of d 3 from volumetric heating calculations without explicit knoweldge of d%. 
This last feature, as well as the similarity of (J2J) and (j7l) is due to fortuitous cancellation; as 
discussed further in section|6j under more general conditions (e.g. higher spatial dimensions), 
this cancellation does not take place and terms containing both d$ and d 3 appear. 
The solution to Equation ^ subject to boundary condition ^ is 



T 



1 4e 



2 5 72 fc 
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— x 



+ d x k - 2d'k 2 



Comparison of this solution to LVDSMC simulations away from the Knudsen layer allows 
us to calculate the coefficient d' 3 . In this work, we extract the value of d' 3 from the slope of 
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T ^ = °^8k' 



di 
2 



(9) 



as a function of k for k — )■ 0. 



4 Computational method 



The Low- variance Deviational Simulation Monte Carlo (LVDSMC) method [T2 | [131 H5 | H"8 | 19] 
efficiently simulates pU] the Boltzmann equation 



written here in the absence of external body forces, by simulating only the deviation f d = 
f — f MB from an equilibrium state / MB . Here, / = f(x, c, t) is the single particle distribution 
function [16] . This approach results in a greatly reduced level of statistical uncertainty for 
low signal problems compared to the standard DSMC [11] approach and is therefore well 
suited to the present application. 

Volumetric heating is modeled by simulating the equation 



dt dx 



df 
dt 



Q (2c 



,.2 



where 



J coll ^0 \6Cq 



f = ^e-H^o (12) 
and po = mn o is a reference mass density. More details on the simulation of the additional 



term on the right hand side can be found in section 

In the versions implemented here, equilibrium is described by a Maxwell-Boltzmann 
distribution 

r B = -^ex P P C -" MB " 2 y (13) 
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based on local (cell-based) mass density p MB , velocity m mb , temperature T MB , and most 
probable velocity c MB = y / 2RT MS . Because f d can take positive and negative values, it is 
represented by signed (or deviational) particles. 

As in the DSMC approach, LVDSMC solves the Boltzmann transport equation through a 
time-splitting approach using a timestep At. The associated advection and collision substeps 
are described below. 



4.1 Advection substep 

During the advection substep, particles move according to the standard DSMC procedure 
(i.e. for particle i, X{(t + At) = Xi(t) + CjAt), with additional particles generated at the 
boundaries and cell interfaces. Each of these additional generation steps are implemented 
by drawing particles from differences of fluxal distributions. 

At a stationary boundary, particles are generated by sampling from 

c ■ n {p B <f) B - f MB ) AAAtd 3 c, (14) 

where AA is the surface area element at the boundary, / MB is the equilibrium distribution 
in the cell adjacent to the boundary, and cf> B is the "boundary distribution" given by 



-c 2 /c 2 



13 



7T 3 / 2 c| 



(15) 



where the c B = a/ 2RT b ; the "boundary density" p B is evaluated from the mass conservation 
statement 

p B f c-n(f) B d 3 c=- f c-nf MB d 3 c. (16) 

J cn>0 J cn<0 
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Particles are also generated at the cell interfaces to account for the spatial discontinuities 
in f MB [T3l [T5] ; they are sampled from 

c • n (f™ B - f™ B ) AA int Atd 3 c, (17) 

where AA mt is the area of the interface, f± B are the equilibrium distributions in adjacent 
cells, and n points from f™ B to f+ B . 



4.2 Collision substep 

The collision substep treatment is based on published LVDSMC implementations [151 E], suit- 
ably modified to include the effect of volumetric heating. We first discuss the BGK collision 
operator and the corresponding volumetric heating implementation; the hard-sphere case 
follows. Due to the small deviations from equilibrium, here we consider the linearized form 
of these collision operators; methods for simulating the corresponding non-linear versions 
can be found in [TBI EH ED] • 



4.2.1 BGK model 

In the case of the BGK model, the collision operator is given by 



df 

at 



f-f 



loc 



coll 



T 



where f loc is the local equilibrium distribution given by 



./ 



loc 



p(x,t) 



[2irRT(x,t)) 



3/2 



exp 



|c — u(x,t)\\' 
2RT{x,t) 



(19) 



where p(x,t), u(x,t) and T(x,t) are the local mass density, flow velocity and temperature. 

Using the approach of Ref. [15], the collision step for the BGK collision operator is 
written as 



df 

at 



At 



At = — [f oc - / MB ] - Af 



At 



coll 



deletion 



(20) 



generation 

where A/ MB is a shift in the equilibrium state. The terms above represent a source term for 
generating new particles, and a sink term for deleting existing particles. It can be shown 
[15] that the generation term is eliminated for linear problems when the equilibrium state 
(for each cell) is shifted according to 



Pmb 
I^mb 
T 



(t + At) 



Pmb 


(tH 


At 


P 


Pmb 






U 




T 

- 1 MB 




T 


T 


— T 

- 1 MB 



(t). 



(21) 



This results in a substantial simplification to step (|20j), which reduces to the very simple 
operation of randomly deleting particles with probability At/r. 
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Because the above method simulates a local equilibrium f MB that is updated in the course 
of the simulation, the heat generation term can be introduced directly (and analytically) into 
the algorithm using 

Q = Po^t (\rTmb ] , (22) 



dt \2 

which results in the following update for the temperature parameter of the equilibrium 
distribution 

2QAi 



AT,, 



3p R 



every timestep. 



(23) 



4.2.2 Hard Sphere model 

The hard sphere collision operator is given by 
'df 



dt 



J coll 



1 

m 




a 



(/ /* - ff*) 



c*||d 3 c* d 2 n, 



(24) 



where primes denote post-collision values and f2 is the spherical angle. The collision step 
for this approach 



\df d ~ 






-s 

coll J_ 



(25) 



v 

generation 



deletion 



is processed as a series of Markov particle generation and deletion steps as proposed by 
Wagner [10]; the specific algorithms employed are discussed in detail in Refs. [TOj E]. In the 
above, 



^ (1) (c,c*) 
u{c) 



0- 2 Pmb 



irmc MB \\c - c 
2 ||c-c,||f 



cxp 



[(c 



71(7 



m 



e-e / 1 



\\c 

L MB I I 



erf(0 



(26) 
(27) 
(28) 



where i = \ \c - u MB ||/c MB . 

In this approach, / MB is not updated during the collision step because the hard-sphere 
simulation algorithm used here is based [9l[T9] on the fixed global equilibrium distribution /°. 
However, to improve accuracy for the low values of Kn considered here Q we have developed 

-'^As shown in |15J . due to the increasing importance of the local equilibrium distribution as Kn — > 0, 
LVDSMC simulations with a variable equilibrium distribution significantly outperform their counterparts 
with fixed equilibrium distribution, because they can be set up to track the local equilibrium distribution 
and thus minimize the number of particles required for the same solution fidelity. 



7 



a special algorithm which uses an equilibrium distribution (/ MB ) that is not updated during 
the collision step but is, however, spatially dependent. In order to determine a suitable form 
of f MB , at the early stages of the simulation, this distribution tracks the local equilibrium 
distribution (similarly to the BGK algorithm described above) using an iterative algorithm 
in which p MB and T MB are taken from the solution at the previous iteration, while the velocity 
u MB is taken to be zero. This process is started with / MB = f° and iterated until / MB no 
longer changes appreciably, which usually takes less than 2 iterations. 

The uniform heat generation is implemented in this case by generating particles from the 
distribution 



or 

dt 



heat 



Q 



2f 
3 eg 



1 / C 



(29) 



Algorithms for efficiently sampling from distributions of the form (29 ) are described elsewhere 

names]. 



5 Results 

Numerical simulations of the uniform heat generation problem were performed in order to 
extract the second-order jump coefficients by comparing the calculated steady centerline 
temperature T(x = 0) with the prediction of equation (|8j at x = 0. 

Figure [I] shows our numerical data for — kd 3 and a linear least squares fit passing through 
the origin based on the data for k < 0.06, and the values d\ = 1.30272 for BGK and 
d\ = 2.4001 for the hard sphere gas [T7j. These fits yield d' 3 = — 1.4 for BGK and —3.1 for 
the hard sphere model; the fit quality demonstrates that the leading order term is indeed k 2 . 
The contribution of higher order terms starts to be noticable as k increases. Incidentally, 
the complementary analysis of the companion paper [7j, based on a finite difference analysis 
of the Knudsen-layer problem of the linearized Boltzmann equation, yields d' 3 = —1.4276 for 
the BGK model and d' 3 = —3.180 for the hard-sphere model. 

Figure [2] shows the temperature field for the hard sphere case with Kn = 0.05 (equivalent 
to k = 0.0443) using the value obtained above (namely d! 3 = —3.1) demonstrating excellent 
agreement everywhere except in the Knudsen layer close to the boundary, as expected. By 
comparing the first- and second-order jump theories, it is clear that the second-order jump 
theory provides an improvement over the existing first-order theory, already at Kn = 0.05. 
For Kn = 0.1 (Figure [3]), the error in the first order solution is quite large, while the second- 
order solution is considerably more accurate, provided that the existence of the Knudsen 
layers for a large part of the domain is accounted for. 



6 Discussion 

Using LVDSMC simulations, we have extracted the second-order temperature jump coeffi- 
cient for a hard-sphere and a BGK gas in the case that the Navier-Stokes-limit behavior 
is captured by an inhomogeneous heat conduction equation, such as the one appearing in 
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k 

Figure 1: Fits used to extract the second-order jump coefficient d' 3 for the hard sphere and 
BGK collision models. 




Figure 2: Second-order temperature jump solution (Equation (|8])) to the uniform heat 
generation problem with Knudsen number Kn = 0.05; simulation results (symbols) are 
compared to the first- (dashed line) and second-order (solid line) jump theories (equation 
(Jsj) with d' 3 = and d' 3 = —3.1, respectively). 
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Figure 3: Second-order temperature jump solution (Equation (|8|) to the uniform heat gen- 
eration problem with Knudsen number Kn = 0.1; simulation results (symbols) are compared 
to the first-order (dashed line) and second-order (solid line) jump theories (equation Q with 
d' 3 = and d' 3 = —3.1, respectively). 
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the presence of constant volumetric heating. Our results have been validated by a com- 
panion paper which provides a deterministic calculation of the same coefficient through a 
rigorous asymptotic analysis of the Boltzmann description of a mathematically equivalent 
problem, namely that of a quescient gas confined between two parallel walls whose tempera- 
ture changes linearly (increases or decreases) in time at a constant (and small) rate. Due to 
the time-dependent nature of the latter problem, the analysis in the companion paper goes 
beyond the asymptotic theory for steady problems [IB] : this also explains why the presently 
calculated jump coefficient (d' 3 ) is not equivalent to the one (d^) obtained by the steady 
asymptotic analysis of Ref . [16] . 

Equation ^ and boundary condition ^ can be generalized to two and three-dimensional 
steady problems as long as the heat generation in the gas is uniform in space and constant 
in time. Specifically, for a quiescent gas, the governing equation and boundary condition in 
this case become 



V 2 T 



5e 



472 k 



(30) 



and 



f\ - T B = (d 1 +d<r ) Kk)k 




V Z T- 



dn 1 



(31) 



respectively. Here R/L is the mean boundary curvature and d^ = 1.82181 for the BGK model 
[16] ; for the hard-sphere gas the value of d§ is unknown. The temperature jump coefficient 
<i 3 for the hard sphere gas, as well as a general second-order slip description of unsteady 
problems remain unknown and will be the subject of future work. 
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